
setwd("C:/Users/Mike/Dropbox/To do/Chapter 1/R")

require(ape)
require(phytools)
require(geiger)
require (mvMORPH)
require(devtools)
require(nlme)
require(qpcR)
require(reshape)
require(ouch)
require(calibrate)
require(surface)


tree<- read.nexus("Chapter 1_trees.nexus")
dat<- read.csv("Chapter 1_data_all PCs.csv", row.names=1)



##Ecology###
eco <- dat[,7]
names(eco) <- dat$Species


TreeOnly <- setdiff(tree$tip.label,rownames(dat)) #checks to make sure the tips of the phylogeny matches the morpho data

DataOnly <- setdiff(rownames(dat), tree$tip.label) #checks to make sure the tips of the phylogeny matches the morpho data

#########################################
#Prune Tree
########################################

phyloTime <- drop.tip(tree,TreeOnly)

phyloTimeLadderized <- ladderize(phyloTime) 
phyloTimeLadderized <- rescale(phyloTimeLadderized, "depth", 1)
plot(phyloTimeLadderized, cex=0.5, no.margin = T) #Plot ladderized and rescaled tree
add.scale.bar() # Add a simple scale bar indicating the scale for the branches in your tree



###model surface + piscivore convergence####
family<- as.vector(dat[,12])  # Make a vector of discreate trait one.
names(family) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(family))],sort(unique(family)))          

mtrees5<-make.simmap(phyloTimeLadderized,family,model="SYM", nsim=1000)

###model for surface####
family<- as.vector(dat[,11])  # Make a vector of discreate trait one.
names(family) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(family))],sort(unique(family)))          

mtrees4<-make.simmap(phyloTimeLadderized,family,model="SYM", nsim=1000)

####model for trophic ecology####
discTraitOne <- as.vector(dat[,6])  # Make a vector of discreate trait one.
names(discTraitOne) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(discTraitOne))],sort(unique(discTraitOne)))          

mtrees<-make.simmap(phyloTimeLadderized,discTraitOne,model="SYM", nsim=1000)


###Model of piscivores###
family<- as.vector(dat[,10])  # Make a vector of discreate trait one.
names(family) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(family))],sort(unique(family)))          

mtrees3<-make.simmap(phyloTimeLadderized,family,model="SYM", nsim=1000)

###Model for families###
family<- as.vector(dat[,8])  # Make a vector of discreate trait one.
names(family) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(family))],sort(unique(family)))          

mtrees1<-make.simmap(phyloTimeLadderized,family,model="SYM", nsim=1000)

####model for continents####
cont<- as.vector(dat[,9])  # Make a vector of discreate trait one.
names(cont) <- rownames(dat) # Name the vector elements.
cols <- setNames(palette()[1:length(unique(cont))],sort(unique(cont)))          

mtrees2<-make.simmap(phyloTimeLadderized,cont,model="SYM", nsim=1000)


 plotSimmap(mtrees[[1]],cols,fsize=0.3,ftype="i")
          add.simmap.legend(colors=cols,prompt=FALSE,x=0.5*par()$usr[1],
          y=30.0*(max(nodeHeights(phyloTimeLadderized))))

pd<-describe.simmap(mtrees,plot=FALSE)
          ## now let's plot a random map, and overlay the posterior probabilities
pdf("ecology_simmap_family.pdf")
plotSimmap(mtrees[[1]],cols,fsize=0.3,ftype="i")
          
nodelabels(pie=pd$ace,piecol=cols,cex=0.5)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.5*par()$usr[1],
y=30.0*(max(nodeHeights(phyloTimeLadderized))))

dev.off()

data<-read.csv("PC1.csv", row.names=1)


OUsurface_pisc<-  mvOU(mtrees5[[1]], dat[], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)

OUsurface<- mvOU(mtrees4[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)

OUtrophic<-mvOU(mtrees[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)


OUpiscivore<- mvOU(mtrees3[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)



OUcont<- mvOU(mtrees2[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)



OUfam<- mvOU(mtrees1[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("OUM", "OU1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)


BM<-mvBM(mtrees2[[1]], dat[,c(1,2,3,4)], error = NULL, model = c("BMM", "BM1"), param = list(sigma = NULL,
          alpha = NULL, vcv = "mvmorph", decomp = c("diagonal")), method = c("rpf", "sparse", "inverse",
          "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B",
          "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL,
          diagnostic = TRUE, echo = TRUE)

EB<- mvEB(mtrees2[[1]], dat[,c(1,2,3,4)], method="pic",  param=list(low=log(10^-5)/10))

results<- list(OUtrophic, OUcont, OUfam, BM, EB, OUpiscivore, OUsurface, OUsurface_pisc)

aicw(results, aicc=TRUE)

